home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / toolbox / besselj.r < prev    next >
Text File  |  1994-09-23  |  8KB  |  301 lines
  1. ////-------------------------------------------------------------------//
  2.  
  3. //  Syntax:  besselj ( ALPHA, X )
  4.  
  5. //  Description:
  6.  
  7. //  Besselj is the bessel function of the first kind of order ALPHA.
  8. //  Either ALPHA or X may a vector in which case a vector result of
  9. //  the same dimension is returned.  If both ALPHA and X are vectors
  10. //  a matrix result of dimension length(ALPHA) x length(X) is returned.
  11. //  If either ALPHA or X is a scalar, the other argument may be a matrix
  12. //  and a matrix result of the same dimension is returned.
  13.  
  14. //-------------------------------------------------------------------//
  15.  
  16. // See Numerical Recipes in C (second edition)
  17. // Currently only integer order supported.
  18.  
  19. static (besseljn, besselj0, besselj1);
  20.  
  21. besselj = function(alph,x) 
  22. {
  23.   global (pi)
  24.  
  25.   // Checking the class will automatically result in an error if
  26.   // the argument doesn't exist.
  27.   if(class(alph) != "num" || class(x) != "num") {
  28.     error("besselj: argument is non-numeric.");
  29.   }
  30.  
  31.   if(min(size(alph)) > 1 && min(size(x)) > 1) {
  32.     // Both arguments are matrices.
  33.     error("besselj: only one argument may be a matrix.");
  34.   }
  35.  
  36.   if(length(alph) == 1) {
  37.     // x is a matrix (or scalar), alph a scalar
  38.     if(alph == 0 || mod(alph,int(alph)) == 0) {
  39.       return besseljn(alph,x);
  40.     else
  41.       error("besselj: fractional orders not yet supported.");
  42.       // return besseljr(alph,x);
  43.     }
  44.   else if(length(x) == 1) {
  45.     // alph is a matrix (or vector), x a scalar
  46.     y = zeros(size(alph));
  47.  
  48.     // Use integer order routines for integer alph
  49.     inte_ind = find(alph == 0 || mod(alph,int(alph)) == 0);
  50.     frac_ind = complement(inte_ind,1:alph.n);
  51.  
  52.     if(inte_ind.n) {
  53.       for(index in inte_ind) {
  54.         y[index] = besseljn(alph[index],x);
  55.       }
  56.     }
  57.     if(frac_ind.n) {
  58.       error("besselj: fractional orders not yet supported.");
  59.       // for(index in frac_ind) {
  60.         // y[index] = besseljr(alph[index],x);
  61.       // }
  62.     }
  63.     return y;
  64.   else
  65.     // alph and x are both vectors
  66.     y = zeros(length(alph),length(x));
  67.  
  68.     // Use integer order routines for integer alph
  69.     inte_ind = find(alph == 0 || mod(alph,int(alph)) == 0);
  70.     frac_ind = complement(inte_ind,1:alph.n);
  71.  
  72.     if(inte_ind.n) {
  73.       for(index in inte_ind) {
  74.         y[index;] = besseljn(alph[index],x);
  75.       }
  76.     }
  77.     if(frac_ind.n) {
  78.       error("besselj: fractional orders not yet supported.");
  79.       // for(index in frac_ind) {
  80.         // y[index;] = besseljr(alph[index],x);
  81.       // }
  82.     }
  83.     return y;
  84.   }}
  85. };
  86.  
  87. // In these functions x can be a vector (or matrix), but n
  88. // must be a scalar.
  89.  
  90. // No argument checking is performed on static functions.
  91.  
  92. besseljn = function ( n, x ) {
  93.  
  94.   local(abs_x,y,index,index_1,index_2,index_3,arg,p_1,p_2,p_3, ...
  95.         two_over_x,m,y_2,Norm,even,limit,lim_ind);
  96.  
  97.   if(n == 0) {
  98.     return besselj0(x);
  99.   else if (n == 1) {
  100.     return besselj1(x);
  101.   else
  102.     // integer order greater than two
  103.     abs_x = abs(x);
  104.     y = zeros(size(x));
  105.     index_1 = find(abs_x > n);
  106.     index_2 = complement(index_1,1:x.n);
  107.     index_3 = find(abs_x == 0);
  108.  
  109.     index_1 = complement(intersection(index_1,index_3),index_1);
  110.     index_2 = complement(intersection(index_2,index_3),index_2);
  111.  
  112.     if(index_1.n) {
  113.       arg = x[index_1];
  114.  
  115.       p_1 = besselj0(arg);
  116.       p_2 = besselj1(arg);
  117.       two_over_x = 2./arg;
  118.       for(index in 1:(n-1)) {
  119.         p_3 = index * (two_over_x .* p_2) - p_1;
  120.         p_1 = p_2;
  121.         p_2 = p_3;
  122.       }
  123.       if(mod(n,2) == 1) {
  124.         index = find(arg < 0.);
  125.         if(index.n) {
  126.           p_3[index] = -p_3[index];
  127.         }
  128.       }
  129.       y[index_1] = p_3;
  130.     }
  131.  
  132.     if(index_2.n) {
  133.       limit=1.e20;
  134.       arg = x[index_2];
  135.       two_over_x = 2./arg;
  136.       // Downward recurrence from arbitrary starting values.
  137.       // Find a starting upper even index (this give 16 digits
  138.       // of accuracy).
  139.       m = 2*int((n + 16*sqrt(n))/2);
  140.       p_1 = zeros(size(arg));
  141.       p_2 = ones(size(arg));
  142.       Norm = zeros(size(arg));
  143.       y_2 = zeros(size(arg));
  144.       even = 0;
  145.       for(index in m:1:-1) {
  146.         p_3 = index * (two_over_x .* p_2) - p_1;
  147.         p_1 = p_2;
  148.         p_2 = p_3;
  149.         if(even) {
  150.           Norm = Norm + p_2;
  151.         }
  152.         lim_ind = find(abs(p_2) > limit);
  153.         if(!isempty(lim_ind)) {
  154.           p_2[lim_ind] = p_2[lim_ind] / limit;
  155.           p_1[lim_ind] = p_1[lim_ind] / limit;
  156.           Norm[lim_ind] = Norm[lim_ind] / limit;
  157.           if(index < n) {
  158.             y_2[lim_ind] = y_2[lim_ind] / limit;
  159.           }
  160.         }
  161.         even = !even;
  162.  
  163.         if(index == n) {
  164.           // Save the unnormalized answer, continue looping to get
  165.           // normalization.
  166.           y_2 = p_1;
  167.         }
  168.       }
  169.       // Only need 1*j0(x) so subtract off p_2 after doubling.
  170.       Norm = 2*Norm - p_2;
  171.       y_2 = y_2 ./ Norm;
  172.       
  173.       if(mod(n,2) == 1) {
  174.         index = find(arg < 0.);
  175.         if(index.n) {
  176.           y_2[index] = -y_2[index];
  177.         }
  178.       }
  179.       y[index_2] = y_2;
  180.     }
  181.  
  182.     if(index_3.n) {
  183.       y[index_3] = zeros(size(index_3));
  184.     }
  185.  
  186.     return y;
  187.   }}
  188. };
  189.  
  190. besselj0 = function(x) {
  191.   local(abs_x,index_1,index_2,y,c_1,c_2,arg_1,arg_2,p_1,p_2);
  192.  
  193.   abs_x = abs(x);
  194.  
  195.   index_1 = find(abs_x < 8);
  196.   index_2 = complement(index_1,1:x.n);
  197.  
  198.   // First evaluate using rational approx. for small arguments.
  199.  
  200.   y = zeros(size(x));
  201.  
  202.   if(index_1.n) {
  203.     c_1 = [ 57568490574, -13362590354, 651619640.7, ...
  204.            -11214424.18, 77392.33017, -184.9052456];
  205.     c_2 = [ 57568490411, 1029532985, 9494680.718, ...
  206.             59272.64853, 267.8532712];
  207.  
  208.     arg_1 = x[index_1] .* x[index_1];
  209.  
  210.     p_1 = c_1[1] + arg_1 .* (c_1[2] + arg_1 .* (c_1[3] + arg_1 .* ...
  211.           (c_1[4] + arg_1 .* (c_1[5] + arg_1 .* c_1[6]))));
  212.     p_2 = c_2[1] + arg_1 .* (c_2[2] + arg_1 .* (c_2[3] + arg_1 .* ...
  213.           (c_2[4] + arg_1 .* (c_2[5] + arg_1))));
  214.     y[index_1] = p_1./p_2;
  215.   }
  216.  
  217.   // Now evaluate for large arguments using approximating form.
  218.  
  219.   if(index_2.n) {
  220.     c_1 = [ -0.1098628627e-2, 0.2734510407e-4, ...
  221.             -0.2073370639e-5, 0.2093887211e-6 ];
  222.     c_2 = [ -0.1562499995e-1, 0.1430488765e-3, ...
  223.             -0.6911147651e-5, 0.7621095161e-6, 0.934935152e-7 ];
  224.  
  225.     arg_1 = 8./x[index_2];
  226.     arg_2 = arg_1 .* arg_1;
  227.  
  228.     p_1 = 1 + arg_2 .* (c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + ...
  229.           arg_2 .* c_1[4])));
  230.     p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* (c_2[4] + ...
  231.           arg_2 .* c_2[5])));
  232.  
  233.     c_1 = abs_x[index_2];
  234.     c_2 = c_1 - pi/4;
  235.     y[index_2] = sqrt(2./(pi*c_1)).*(p_1.*cos(c_2) - arg_1.*p_2.*sin(c_2));
  236.   }
  237.  
  238.   return y;
  239. };
  240.  
  241.  
  242. besselj1 = function(x) {
  243.   local(abs_x,index_1,index_2,y,c_1,c_2,arg_1,arg_2,p_1,p_2);
  244.  
  245.   abs_x = abs(x);
  246.  
  247.   index_1 = find(abs_x < 8);
  248.   index_2 = complement(index_1,1:x.n);
  249.  
  250.   // First evaluate using rational approx. for small arguments.
  251.  
  252.   y = zeros(size(x));
  253.  
  254.   if(index_1.n) {
  255.     c_1 = [ 72362614232, -7895059235, 242396853.1, ...
  256.             -2972611.439, 15704.48260, -30.16036606 ]; 
  257.            
  258.     c_2 = [ 144725228442, 2300535178, 18583304.74, ...
  259.             99447.43394, 376.9991397];
  260.  
  261.     
  262.     arg_1 = x[index_1];
  263.     arg_2 = arg_1 .* arg_1;
  264.  
  265.     p_1 = arg_1.*(c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + arg_2 .* ...
  266.           (c_1[4] + arg_2 .* (c_1[5] + arg_2 .* c_1[6])))));
  267.     p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* ...
  268.           (c_2[4] + arg_2 .* (c_2[5] + arg_2))));
  269.     y[index_1] = p_1./p_2;
  270.   }
  271.  
  272.   // Now evaluate for large arguments using approximating form.
  273.  
  274.   if(index_2.n) {
  275.     c_1 = [ 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, ...
  276.             -0.240337019e-6 ];
  277.     c_2 = [ 0.04687499995, -0.2002690873e-3, 0.8449199096e-5, ...
  278.             -0.88228987e-6, 0.105787412e-6 ];
  279.  
  280.     arg_1 = 8./x[index_2];
  281.     arg_2 = arg_1 .* arg_1;
  282.  
  283.     p_1 = 1 + arg_2 .* (c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + ...
  284.           arg_2 .* c_1[4])));
  285.     p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* (c_2[4] + ...
  286.           arg_2 .* c_2[5])));
  287.  
  288.     c_1 = abs_x[index_2];
  289.     c_2 = c_1 - 3*pi/4;
  290.     arg_2 = sqrt(2./(pi*c_1)).*(p_1.*cos(c_2) - arg_1.*p_2.*sin(c_2));
  291.     index_1 = find(arg_1 < 0);
  292.     if(!isempty(index_1)) {
  293.       arg_2[index_1] = -arg_2[index_1];
  294.     }
  295.  
  296.     y[index_2] = arg_2;
  297.   }
  298.  
  299.   return y;
  300. };
  301.